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Abstract. We consider the realization of universal quantum computation through 
braiding of Majorana fermions supplemented by unprotected preparation of noisy 
ancillae. It has been shown by Bravyi [Phys. Rev. A 73, 042313 (2006)] that under the 
assumption of perfect braiding operations, universal quantum computation is possible 
if the noise rate on a particular 4-fermion ancilla is below 40%. We show that beyond a 
noise rate of 89% on this ancilla the quantum computation can be efficiently simulated 
classically: we explicitly show that the noisy ancilla is a convex mixture of Gaussian 
■^j- | fermionic states in this region, while for noise rates below 53% we prove that the 

C*") ■ state is not a mixture of Gaussian states. These results were obtained by generalizing 

concepts in entanglement theory to the setting of Gaussian states and their convex 
mixtures. In particular we develop a complete set of criteria, namely the existence of 
a Gaussian-symmetric extension, which determine whether a state is a convex mixture 



OO 
O. 

, of Gaussian states. 



1. Introduction 

One interesting route towards universal quantum computation is through the realization 
of Majorana fermion qubits. Such Majorana fermion qubits, encoded in pairs 
of nonlocal fermionic zero-energy modes, are believed to be present in various 
quantum systems such as a v = 5/2 fractional Quantum Hall system, p x + ip y 
superconductors, and recently proposed topological insulator /superconductor and 
semiconducting nanowires/superconductor structures. See p] for a review. Braiding 
operations on such Majorana fermion qubits can implement certain topologically- 
protected gates, namely single-qubit Clifford gates [2]. A system of Majorana fermions 
and braiding operations (which are a subset of non-interacting fermion operations) can 
be classically efficiently simulated if it is not supplemented by additional resources. 
Universal quantum computation can be achieved if this supplement consists, for 
example, of either (i) gates which use a quartic interaction between Majorana fermions or 
(ii) a quartic parity measurement or (iii) two ancillae |a 4 ) and \a 8 ) involving respectively 
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4 and 8 Majorana fermions. The advantage of the last realization is that even when these 
ancillae are noisy, one can purify them using braiding operations. This scheme is similar 
to the magic-state-distillation scheme by which Clifford group operations are used to 
distill almost noise- free single qubit 7r/8 ancillae from noisy ones [3] . The noise threshold 
for such distillation schemes, which assume that the distilling gates and operations are 
noise-free is of the order of tens of percents what makes them attractive. But one can 
also ask the converse question: how noisy are these ancillae allowed to be before one 
can efficiently classically simulate the entire quantum computation? This question has 
been addressed in the case of noisy single-qubit magic states and noisy single-qubit gates 

In this paper we consider a similar question for computation using Majorana 
fermions. It is not hard to show (see Section [5]) that if the noisy ancillae are convex 
mixtures of Gaussian fermionic states and the computation involves only non-interacting 
fermionic operations that one can still classically efficiently simulate such quantum 
computation. Hence we set out to develop a criterion that determines whether a state 
is a convex mixture of Gaussian fermionic states in Section [31 We find such criterion 
in the form of a hierarchy of semidefinite programs, similar as for separable states [7]. 
Unfortunately, the computational effort for implementing this general criterion is too 
large to give decisive information for our problem of interest and we have recourse to 
an analytical approach for the noisy \a$) ancilla in Section HI Even though our general 
criterion is not immediately useful for the problem at hand, its generality and similarity 
with separability criteria makes it interesting in itself. 

Our work is different from previous work on criteria which determine whether 
a fermionic state with a fixed number of fermions has a single Slater determinant, 
or whether a fermionic state can be written as a convex combination of states with 
single Slater determinant [U |9]. The important distinction is that we fix only the 
parity of the fermionic state and not the number of fermions. Our goal is to extend 
the class of Gaussian fermionic states in a natural way, by considering states which 
can be written as convex combinations of Gaussian states. A state which is such a 
convex combination of Gaussian states we call convex-Gaussian or 'having a Gaussian 
decomposition'. Even though physical states have a fixed number of fermions, Gaussian 
fermionic states are important approximations to physically non-trivial states, such 
as the superconducting BCS state. Our criterion thus intends to separate Gaussian 
states and convex combinations thereof, from fermionic states with a richer structure 
which cannot be simply viewed as states in which fermions are paired fl Note that the 
question of having a Gaussian decomposition is also different from the question of pairing 
which is analyzed in [TTJ: Gaussian fermionic states can be paired while single Slater 
determinant states cannot. We hope that our results in separating convex- Gaussian 
states from fermionic states with a richer structure may lead to tools for understanding 

J Any Gaussian fermionic state with an even number of fermions can be brought, by fermion- number 
preserving quadratic interactions, to a normal form which is a superposition of states with fixed pairs 
of fermions created from the vacuum, see e.g. [TO] . 
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groundstates of interacting fermion systems beyond mean-field theory. 
2. Preliminaries: Fermionic Gaussian states 

We consider a system of m fermionic modes, with corresponding creation (a\.) and 
annihilation (a*,) operators (A; = 1, . . . m), respecting the Fermi-Dirac anti-commutation 
rules, {aj,ak} = and {aj,a\.} = Sj^I. For systems in which only fermionic parity is 
conserved, it is more convenient to use the 2m Majorana fermion operators defined as 
C2/C-1 = Gfc + % and C2k = ~K a l ~ a k)- The operators are hermitian, traceless 

and form a Clifford algebra C 2m with {cj, c^} = 2^/0. 

Any Hermitian operator X G C 2m that can be written as the linear combination of 
products of an even number of Majorana operators is called an even operator, i. e. 

m 

X = apl + ^ ] i ^ ] 0i a 1 ,a2,...,a 2k c a 1 Ca2 ■ ■ ■ c a 2k i (1) 

k=l I<ai<a2<...<a2fc<2m 

where the coefficients «o and a ftl(0a) ... |aaJi are real, for (c ai c a2 . . . c a2k Y = 
(—l) k c ai c a2 . . -c a2k . The parity of the number of fermions is conserved by the action 
of an even operator X as X commutes with the fermionic number-parity operator 
Call = i m cic 2 . . . c 2m = rGLi(^ ~ 2ajaj). Thus the projector onto a pure state \ip) 
with fixed parity is an even fermionic operator, while has eigenvalues C a n = ±1 
depending on whether the parity of the number of fermions in is odd or even. 

Given an even state g G C 2m , the correlation matrix M is a 2m x 2m real, anti- 
symmetric matrix with elements 

M ab := ^Ti(g[c a ,c b ]), with a, b = 1, . . . , 2m. (2) 

Real, anti-symmetric matrix such as M can be brought to block-diagonal form by a real 
orthogonal transformation R G 5*0 (2m). i.e. 

Let us now define fermionic Gaussian states. Fermionic Gaussian states g are 
even states of the form g = K exp(— i£\ , } . flijCiCj) with real anti-symmetric matrix /3^- 
and normalization K. Hence Gaussian fermionic states are ground-states and thermal 
states of non-interacting fermion systems. One can block-diagonalize /3 and re-express 
q in standard form as 

^ m 

8 = o"^ II Xk ' Where Xk = 1 + lX kC2k-lC2k, (4) 
k=l 

where c = R T c with R block-diagonalizing the matrix fy. The coefficients Aj (which 
can be expressed in terms of the eigenvalues of will lay in the interval [—1,1]. For 

§ A system of m qubits is isomorphic with a system with m fermions and the unitary Jordan- Wigner 
transformation maps each Majorana fermion operator onto a nonlocal product of Pauli operators acting 
on m qubits. 
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Gaussian pure states, Xj G { — 1, 1} so that M T M = I, while for Gaussian mixed state 
M T M < I. In the theory of Gaussian fermionic states, a special role is played by those 
unitary transformations which map Gaussian states onto Gaussian states. These are the 
transformations generated by Hamiltonians of non-interacting fermionic systems, i.e. 
Hamiltonians which are quadratic in Majorana fermion operators. In this paper we will 
refer to these unitary transformations as fermionic linear optics (FLO) transformations, 
as they, similarly as for bosonic linear optics transformations, have the property that 



where U is a FLO transformation and R G SO (2m). Hence the transformations which 
block-diagonalize M (or 0) are FLO transformations. Note that C a n is invariant under 
any unitary FLO transformation U, as UC a \\ = C a \\U. 

It is clear from the standard form, Eq. (]!]), that a fermionic state is fully determined 
by its correlation matrix M. The expectation of higher order correlators for a Gaussian 
fermionic state are efficiently obtained by virtue of Wick's theorem 



where M| aii _ a2p is the sub-matrix of M which contains only the elements Mj k with 
j, k — ai, . . . , a 2p and Pf(.) is the Pfaffian [jj. 

Their efficient description (j2J) combined with the possibility to efficiently evaluate 
the expectation value of observables fl6]) makes fermionic Gaussian states a valuable 
tool for approximating ground-states, thermal states or dynamically-generated states 
using (generalized) Hartree-Fock methods of interacting fermion systems, see e.g. 
[TUl H"2] . Their concise representation together with Eq. (jSJ) also allows for an efficient 
classical simulation of quantum computations that employ only FLO operations and are 
initialized with Gaussian states [131 El E3 US]- 

3. Beyond Gaussian states: Gaussian decompositions 

Given all the applications of Gaussian states, it is natural to try to enlarge this 
set of states to form a convex set, in analogy with extending product states to 
separable states. A first observation in this direction is that mixed Gaussian states 
are straightforward convex mixtures of pure Gaussian states. This is easily seen using 
the standard form fl3J), i.e., any Gaussian state can be written as g = Y[T=i Xk/2 m with 
Xk = I + i(2p k -l)c 2 k-iC2k = Pk(I + ic 2 k-iC2k) + (1 -Pk)(I -ic 2 fc-ic 2 fc) where <p k < 1. 
Note that we get the completely mixed state I /2 m by setting p k = 1/2. The converse, 
however, is not true: the convex mixture of two Gaussian states is in general not a 
Gaussian state [T71 [18] . Let us define 

|| Let A = ( a i,j) be a 2m x 2m anti-symmetric matrix, then Pf(^4) = 
2^7S! S we s 2m sgn(7r) III™, i a 7r(2i-i),7r(2<)) where S 2m is the set of all permutation of 2m symbols. 
The Pfaffian is non-zero only for anti-symmetric matrix of even dimension. 




(5) 



j 




(6) 
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Definition 1 (Convex-Gaussian). An even density matrix g G C 2m is convex-Gaussian 
iff it can be written as a convex combination of pure Gaussian states in C 2m , ^.e. 
g = J^i.Pi^i where o~i are pure Gaussian states, pi > and YliPi = 1- 

Remark 1: Since every even state in C 2m is defined by 2 2m — 1 real coefficients, 
as is clear from Eq. (CQ) fixing the normalization «o = l/2 m , Caratheodory's theorem 
guarantees us that for any convex- Gaussian state there exists a Gaussian decomposition 
with at most 2 2m pure Gaussian states. 

Remark 2: There is both similarity and difference between the set of separable 
states and the set of convex-Gaussian states. For example, for separable states local 
unitary transformations captured by 0(m) parameters relate the extreme points of the 
set, whereas for convex-Gaussian states the extreme points are the pure Gaussian states 
which are related by FLO transformations, captured by 0(m 2 ) parameters. Note that 
the pure states o\ in decomposition of a convex-Gaussian state will typically have a 
correlation matrix M{<Ji) which is block-diagonal in a different {q} basis. If we were to 
restrict ourselves to convex combinations of Gaussian states which are in standard form 
using one and the same set {q}, we would obtain a convex set isomorphic to the set of 
separable states diagonal in the classical bit-string basis. 

We will now prove that for systems of 1, 2 or 3 fermionic modes, all pure even 
states are Gaussian (a different proof of this result can be found in [19]). From this 
observation it is then immediate that any even state g G C 2m , with m = 1,2, 3, is 
convex-Gaussian. In order to prove this result we start with a useful Lemma which 
shows that any fermionic even density matrix has a correlation matrix with eigenvalues 
in the complex [— i, i] interval. 

Lemma 1. The correlation matrix M of any even density matrix g G C 2m has 
eigenvalues ±iAfc with \ k G [—1,1] and k G {l,...,m}. M T M = I iff g is a Gaussian 
pure state. 

Proof. In order to prove this, we use an approach introduced in [20]. Let g G C 2m be 
any density matrix and let {cj} 2 !!\ be the set of Majorana fermions in which M is block- 
diagonal (if the choice for {c{\ is not unique, pick one). Define the FLO transformation 
U k , k = 1, . . . , m as 

U k c 2k Ul = -c 2k , U k c 2k -iUl = -c 2k -i, U k CiUl = Ci. (7) 

Note that the FLO transformation U k induces an orthogonal transformation R k that 
leaves the correlation matrix M of g invariant. Let g k = ^(Qk-i+U k Qk-iUl), with £ = Q- 
The iteration in k has the effect of dephasing g in the eigenbasis of each \C2 k -\C2 k - 
Thus g m is a density matrix which contains only the mutually commuting operators 
ic 2 fc_ic 2 £;, and hence its eigendecomposition involves the eigenstates of these operators 
which are Gaussian pure states. We can thus represent the state g m = ^2 x p x \x){x\, 
where x G {0, l} m and p x > with Y2 x Px = 1- ^ ne correlation matrix M of g m (and 
thus of g) is that of a convex-Gaussian state with pure Gaussian states in the same 
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standard form which implies that the correlation matrix M is that of a Gaussian mixed 
state. Note that the dephasing procedure can only increase the entropy of g. Hence if 
the final state has the property that M T M = I, i.e. corresponds to a pure Gaussian 
state then g itself should have been a pure Gaussian state. Thus M T M = I iff g 

is a pure Gaussian state. □ 

Proposition 1. Any even pure state \ip)(ip\ G C 2m form = 1,2,3 is Gaussian. 

Proof. The statement is trivial for m = 1. Consider m = 2 so that the standard 
form of any even pure state can written as some linear combination of /, quadratic 
terms ic 2 / c _ic 2 / c and C a n. The dephasing procedure in the proof of Lemma [1] leaves 
any such even density matrix invariant - remember that the operator C a \\ is invariant 
under FLO transformations. As the output is always convex-Gaussian, then the input 
state was a pure Gaussian state. The proof for m = 3 follows the same structure 
as for m = 2, but employs the additional observation that C a \\\il)){ijj\ = ±\ip)(ip\ for 
even pure states. The standard form of an even density matrix G C 6 can be 

expressed as = al + (3C & ii+iJ2k7kC2k-iC2k + J2i<j<k<iVijkiCiCjC k ci. The condition 

that C a \\\i/;) = ±\i/j)(i/;\ implies that the quartic terms only have paired correlators 
c 2 k-iC2k'- this is because C a n interchanges the quartic and quadratic terms. Thus \i>){i^\ 
has quartic terms which only involve commuting operators c 2 fc_ic 2 fc. As argued in the 
proof of Lemma [H \if))(ip\ is then a convex mixture of pure Gaussian states, and the 
result follows. □ 



As for separable states j2Tj (22J [23] , the volume of convex-Gaussian states is different 
from zero for any m: any state g G Ci m is convex- Gaussian if it is "close enough" to 
the maximally mixed state J/2 m , which is Gaussian. The following Lemma proves this, 
using arguments similar as those in Ref. 



Theorem 1. For every even state g G €2™, there exists a sufficiently small e > such 
that g e = eg + (1 — e)I/2 m is convex-Gaussian. 

Proof. Consider the set {<t(A,7t)} of pure Gaussian states, with each a(X, tt) G C 2m 
defined as 

7T) = — Y[(I + i\ k C^(2k-l)C^2k)), (8) 
fc=l 

where A G { — l,l} m and n G S£ m , with the set of permutations with positive 
signature (these correspond to FLO operations). There are (2m)\/m\ pure states in this 
set and we will show that they form an over-complete basis for C 2m . We prove this by 
showing that any Hermitian g can be written as 

(X, P i) 
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where w,x n > and c > 0. Furthermore, we can write 

(A,7T) — — ' 

A,7T 

If we find the decomposition in Eq. fl9]) with a certain c (which one could try to minimize 
in our construction of the decomposition), it then follows immediately that j^Q+ ^ ^ 
is convex-Gaussian. Taking e = ^ we obtain our result. One could in principle upper 
bound c as a function of m. Eq. fl9]) can be proved by considering how to express a single 
correlator a ai ^.. )a2k c ai c a2 . . . c a2k in terms of Gaussian states and /. Let us call the subset 
{ai, . . . , a 2 k} = S and the coefficient a ai) ... )02fc = a. Let C5 = c ai c a2 . . . c a2k . We will use 
Gaussian mixed states £(A,7r) = FlfcLiC + i^kCn(2k-i)Cir(2k)) where A, G [-1, 1]; such 
states can in turn be expressed as convex mixtures of pure states <r(A, 7r) for Aj G {—1, 1} 
as in Eq. (jHJ). We have 

i fc TrCs£(A, 7r) = n fc:7r (2fc_i) i7r (2fc)G5Afc- (10) 

By choosing X k G {—1, 1,0} we can make sure that (i) |a|£(A,7r) = w^ w ^(X,tt) has the 
correct expectation a for the correlator C,s (note that the sign of A& can be chosen to fit 
the sign of a) and (2) £(A, it) has no higher-order correlators (by choosing some Afe = 0). 
The state £(A,7r) will then also have lower- order non-zero correlators Cs> for S' C 5 1 , 
including 5" = corresponding to /. Thus we repeat the procedure with these various 
lower-order correlators, i.e. we represent them as a Gaussian mixed state £(A, 7r) plus 
additional correlators corresponding to subsets of S' . Proceeding this way, we remove 
all correlators except for / and end up expressing aCs = J2\n w (\ 7r)£(^> 77 ) — ^ ^ or 
p = 2~ m J2\n w ^- This procedure is probably far from optimal, i.e. gives rise to a 
large c. One can optimize this by applying the procedure to a density matrix g from 
which we repeatedly take out Gaussian mixed states, i.e. g — >■ g — which remove the 
current highest-order correlator and the Gaussian state £ is chosen to be the one with 
minimum (but non-negative) weight Trw£. □ 

Now that we have established the essential features of the set of convex- Gaussian 
states, the natural question that follows is how to determine whether a state has a 
decomposition in terms of Gaussian pure states. As for separability, there are two 
sides to this question. One is to find a Gaussian decomposition, a problem which 
is likely, as in the entanglement case, to be computationally hard (NP-hard in the 
dimension of the space iV = 2 m ) in general (see Proposition [2]). The reverse question is 
to find a criterion which establishes that a state is not convex-Gaussian. Note that 
the goal is to find a criterion which acts linearly on the pure Gaussian states, so 
that it can naturally be extended to convex mixtures thereof. The hermitian operator 
'=1 Ci®c.i G Cim ®C>im-, introduced in Ref. [14J, will be useful in this context as it 
has been shown to lead to a necessary and sufficient criteria for a state to be Gaussian: 

Lemma 2 (Bravyi [H]). An even state g G C 2m is Gaussian iff [A, g ® g] = 0. 
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The operator A has the important property of being invariant under U ®U where U is 
any FLO transformation. The action of A on g <g> g can be appreciated by considering 
the trace norm ||.||i of the positive operator Ag ® gA: 

HA^AHi = 2m-J2(Tric a c b gf = 2m-^(M ab f = 2m—TrM T M, (11) 

a^b a,b 

where M is the correlation matrix of the state g. For a pure Gaussian state g = a, we 
have M T M = I implying that A a® a = 0. For a mixed Gaussian state or non-Gaussian 
state M T M < I, see Lemma [U and hence ||A^>® gA\\i > 0. These observations are 
collected in the following corollary. 

Corollary 1. For an even state g G C2 m , Ag®g = 0iffgisa pure Gaussian state. 

Corollary [T] says that the state of two copies of pure Gaussian states is contained in 
the null-space of the operator A. In Lemma [4] in the Appendix A we will prove that the 
null-space of the operator A is contained in the symmetric subspace Sym 2 (C 2m ) spanned 
by vectors which are invariant under the SWAP operator P, i.e. = P\i(j). 

In particular, the null-space of A is spanned by states \ip,ip) where ip is Gaussian, 
while the symmetric subspace is spanned by \<f>, (ft) for any 0. We will call the null- 
space of A the Gaussian-symmetric subspace. As A is a sum of mutually commuting 
hermitian operators q £g> q with /i« = ±1 eigenvalues, the projectors on the eigenstates 
are P^ = ^^11^(1 + fiiCi <8> q) with eigenvalue Y^i^i °f A. Hence the Gaussian- 
symmetric subspace is spanned by ( 2 ^) orthogonal eigenvectors P^ with the property 
that fii = 0. Clearly, the dimension of the Gaussian-symmetric subspace ( 2 ™) is 
smaller than the dimension of the symmetric subspace Sym 2 (C 2m ) which is ( 2 2 +1 ) for 
any m > 1. 

In [7J [21] the authors established a series of tests for separability based on the 
existence of a symmetric extension of any separable state. The usefulness of these tests 
rely on the fact that they correspond to semi-definite programs which for small numbers 
of extensions can be implemented numerically. Furthermore, these tests are complete, 
in the sense that an entangled state does not have a symmetric extension to an arbitrary 
number of parties in the extension. Here we formulate a similar series of extension tests 
which are all passed by convex-Gaussian states, but non-convex-Gaussian states fail in 
some of them. Our criterion is based on the following simple observation: Let a density 
matrix g G Cim be convex-Gaussian, i.e. we can write g = ^2iPiO~i with o~i pure Gaussian 
states, then there exists a symmetric extension g ext G Cf™, namely £ ext = J2iPi a ? n i 
which is annihilated by A acting on any pair of tensor-factors, and Tr 2v .. n -ig ex t — Q- 
This immediately leads to the following feasibility semi-definite program: 
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Program 1. Input: g G C 2 ra o,nd an integer n > 2 
Body: Is there a g ext G Cf™ 
s-t. Tr 2 ,...,n6ext = Q 

Tr g ext = 1 
A fe 'W = 0,VM/ 

£ext > 

Output: yes (provide g ext ) or no 
Note that we do not need to enforce any permutation-symmetry on the extension 
as the condition that A k ' l g ext = where A k)l = J2A c i)k ® {<k)h acting on tensor-factors 
k I, forces, by LemmaHJin Appendix A, the extension to be in the symmetric subspace 
of all parties. It is easy to see that any state g has an extension supported in the 
symmetic subspace, hence the constraint imposed by A is crucial. 

We say that a state g has a Gaussian-symmetric extension if there exists a n for 
which the Program [1] returns yes. An explicit construction of the SDP for the case 
n = 2 is given in Appendix B, and some numerical results are discussed in the following 
section. Here we will prove that only convex- Gaussian states have a Gaussian-symmetric 
extension to an arbitrary number of spaces, thus showing that these criteria form a 
complete hierarchy for deciding membership in the set of convex- Gaussian states. 

Theorem 2. An even state g G C 2m has a Gaussian-symmetric extension to an arbitrary 
number of parties if and only if g is convex-Gaussian. 

Proof. One direction is immediate, and has been already alluded to in the formulation of 
the SDP. Assume that g is convex-Gaussian, i.e. g = J2iPi a i where ex, are pure Gaussian 
states and {p^ is a probability distribution. Then the extension g cxt = J2iPi a f n nas 
the property that A k,l g™ xt = where A fc '' acts on the tensor-factors k and / by Corollary 

m 

To prove the other direction, we will invoke the quantum de Finetti theorem |25j . 
Let g™ xt > be the extension of g to Cf™ such that A k ' l g r ^ xt = and Tr 2i ... jn £?e X t = 6- 
Using Lemma H] this shows that g™ xt is in the symmetric subspace. Let g™ 2 be the 
reduced density matrix for the first two tensorfactors, i.e. g" 2 — Tr 3v . ., n Qext- Then 
Theorem II. 8 in |25j tells us that 

1 1 0i,a - ^2la(n)r a (n) ® r a (n)||i < e m (n), (12) 

a 

with e TO (n) = 4x ^ . Here r a (n) G C2m are density matrices and j a (n) are probabilities 
summing up to 1. We have 

1 1 A la(n)T a (n) ® r a (n)A \\i = \ \ A(^ 2 - ^ la(n)T a {n) <g> r a (n))A | |i 

a a 

< ||A||ie m (n) = e m (n), (13) 

with ||A||f = 4(m + l) 2 ((^ 1 )) 2 - For fixed m, in the limit n — > oo the upper-bounds 
e m (n),e m (n) -> 0, and by Eq. (JT2J) g\ 2 -)> Yl, a la{n)T a {n) (8) r a (n). In this limit, Eq.(H3D 

% Adding the permutation-symmetry would of course leave the program in SDP form. 
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implies, using the observation in Eq.( lTT|) . that for each a either j a (n) — > or T a {n) 
converges to a Gaussian pure state. Hence for n — )■ oo the state g\ 2 converges to a 
convex mixture of two copies of Gaussian pure states, and its reduced density matrix 
converges to a convex-Gaussian state. □ 

Remark 3: More direct proofs (by establishing a Gaussian-fermionic de Finetti 
theorem directly using the existence of an extension in the Gaussian-fermionic subspace) 
and quantitative versions of this result should be possible. 

Note that convex-Gaussian states always have an symmetric extension which is 
separable in relation to all tensor-factors. This property can be used to provide an 
alternative formulation of the necessary and sufficient criterion for a state to be convex- 
Gaussian: 

Proposition 2. A state g G C 2m is convex- Gaussian iff there exists an extension 
g ext G C 2m <g) C 2m that is a feasible solution to Program U\ -input g and n = 2- with 
the additional constraint that g ext is separable between the tens or- factors. 

Proof. As before one direction is immediate. For the other direction, assume that there 
exists a state g ey ± G C 2m ®C<i m which is separable, i.e. it is of the form g cxt = YliiPi T i®^ii 
with {rj} and {^} sets of density matrices in C 2m and {p^} a probability distribution. 
From the constraint Ag ext = 0, it follows that g is in the symmetric subspace, t% = q, 
and furthermore Tj must be a Gaussian pure state. Therefore, g = Tr 2 ^iPi r i ® r « is 
convex-Gaussian. □ 

Including this extra separability constraint in Program [1] breaks its SDP structure 
- checking separability can in itself be cast as a SDP [21] , but this would give a nested 
SDP structure to the program deciding if a state is convex-Gaussian. A direct semi- 
definite relaxation of this criteria is obtained by employing the partial-transposition test 
for separability [26, 127] , leading to the following SDP: 

Program 2. Input: g G C 2m 

Body: Is there g ext G C 2m <8> C 2m 
s.t. Tr 2 g e xt = g 
Tig cxt = 1 

J^Qext = 
£ext > 

eSt > 

Output: yes (provide g ext ) or no 

Here T\ is the transposition map on the first tensor-factor. As for the previous program, 
a negative answer means that the state is not convex-Gaussian. For the converse, it is 
likely that the PPT criterion together with the constraint that the extension lives in the 
Gaussian-symmetric subspace is not sufficient to enforce separability (e. g. there exists 
symmetric bound-entangled states [28]), hence a positive answer is non-decisive, but we 
leave this as an open question. 
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4. Applications: depolarized |a8) 

As mentioned in the introduction, one way to turn v = 5/2 topological quantum 
computation into an universal quantum computation scheme, is to assume that one 
has access to the auxiliary states (04) G C4 and |a§) G C$, as shown in [2]. In the present 
section we want to assess the amount of noise that can be added to these extra resources 
before they turn convex-Gaussian. From Proposition [1] we know that any pure state in 
C4 is Gaussian, and therefore any noisy version of ^4) is convex-Gaussian. We thus 
concentrate on the state |ag) when undergoing depolarization. We consider the state 

Q aS (p) = (l-p)\a 8 )(a 8 \+pI/16, (14) 

with < p < 1 and 

1 



|«8)(«s| 



16 



{I + S 1 )(I + S 2 )(I + S 3 )(I + Q), 



(15) 



with S 1 = -C!C2C 5 c 6 , #2 = -c 2 c 3 c 6 c 7 , 5*3 = -C1C2C3C4 and Q = c^c^c^CqC^. Ref. [2] 
has shown that g aa (p) can be distilled for p < 40% (corresponding to e 8 = 0.38). We 
would like to determine the noise threshold p* such that for g a %{p > p*) is convex- 
Gaussian. Note that g as is Gaussian only for p = 1. The results in this Section are 
summarized in Fig. [TJ 



p=0.40 p=8/15 



Classical Simulation 
Convex Gaussian 




p=8/9 



Figure 1. Properties of the state g as (p) = (1 — p)\ a 8)( a &\ +p//16. Our results show 
that p t , the noise threshold above which g a g(p) is convex-Gaussian is in the interval 
[ji;, |], depicted in grey. In this interval we do not yet know whether g as is convex- 
Gaussian. For p < 0.40, Ref. [2] shows that g as can be distilled by noiseless braiding 
operations to enable quantum computation with Clifford group elements (for quantum 
universality one also needs a sufficiently low-noise |a4) ancilla). Our results from the 
SDP on the existence of a Gaussian-symmetric extension of g as (p) f° r P ^ 0.25 do not 
provide additional information beyond these analytical results. 



In order to give an lower-bound on p* we implemented the SDP Program [T]for fixed 
n = 2 (see Appendix B) with the aid of the cvxopt library for Sage [29J. The numerical 
results indicate that > 0.25, as it is always possible to find a Gaussian-symmetric 
extension for Qasip) with n = 2 for noise rates above this value. Our implementation of 
the feasibility SDP takes approximately one day for each value of p in a computer with 
2.3GHz processor, and takes about 30GB of RAM memory. The implementation of the 
SDP with an additional PPT condition on the extension is likely to provide stronger 
results, but it is more computationally intensive. Given these limitations, we provide an 
alternative analysis of the properties of g a8 (p) using the notion of witnesses: hermitian 



The Power of Noisy Fermionic Quantum Computation 



12 



operators which have expectation values in a restricted range (e. g. non-negative for all 
Gaussian states) whereas expectation values on non-Gaussian states can extend beyond 
this range (e. g. be negative). The existence of such witness operators is guaranteed 
by the fact that they represent separating hyperplanes separating the convex- Gaussian 
states from all even hermitian operators in €2™, (in addition, similar as in [21], witnesses 
could be constructed based on a negative output of feasibility SDP). 

Clearly, such hermitian witness operators, say W, should include terms which are 
quartic or higher-weight correlators, as the expectation value of quadratic Hamiltonians 
only depends on the correlation matrix of a state. For £> a8 (p) the obvious choice for such 
a witness operator is W = |a 8 )(a 8 |, for which TrWg a s(p) = 1 — Let us understand 
how we can bound min^ g TTW\ip g )(ip g \ = \(ip g \ag)\ 2 where ip g is any Gaussian pure state. 
We will use the fact that the correlation matrix of |a 8 ) is the null matrix, i.e. for all 
Majorana operators q,Cj, Tr CiCj\a%)(a%\ = 0, which follows directly from Eq. (ITS"]) . We 
will now prove that \(ip g \a 8 )\ 2 < ~ which shows that for p < ^ ps 0.53, the state Q a %{j>) 
is not convex-Gaussian. The following proposition proves our claimed bound and can 
be intuitively understood by interpreting the state |a 8 ) as a maximally entangled state 
while Gaussian states are product states: 

Proposition 3. For all Gaussian pure states \ip g ), we have \ (ijj g \as) | 2 < \- 

Proof. Let c\ and C2, C\C2 = — C2C1, be two arbitrary Majorana fermion operators, from a 
complete set {c 2 k-i, c 2 k}i=i (note the 5*0(8) freedom in choosing these). Let \x) denote 
an eigenvector of the mutually commuting operators iciC2, \cik-\Cik for k = 2, ... 4 where 
the bits (0 or 1) of x label the eigenvalues (resp. +1 and —1) of these operators. The 
Gaussian states \x) form a basis for the 4-fermion (or 4-qubit) space, hence we can write 

K) = 5Le{o,i} 4 a x\ x ) witn J2 X \ a *\ 2 = !■ Then 

Tricic 2 |a 8 )(a8| = \ a 0y\ 2 ~ l" 1 ^ 2 = °" ( 16 ) 

j/e{o,i} 3 2/e{o,i} 3 

This, together with normalization, implies that | «o y | 2 = J2 y \ a iy\ 2 = |' hence for all 
y, \ctoy\ 2 and \cti y \ 2 are less than or equal to 1/2. □ 

A tighter bound than the one in Proposition[3]cannot be excluded if one uses further 
properties of the state |a 8 ) - such tighter bound would lead to a reduction of the grey 
area in Fig. [U 

To give an upper-bound on p* we construct explicit Gaussian decompositions of 
£> a8 (p) for large values of p. To do this we consider a subset of the convex- Gaussian 
states that share some key properties with Qasip), namely: (i) zero correlation matrix, 
and (ii) only a small fraction of all possible correlators has non-zero coefficients. 

To exploit these symmetries we define three types of states: 

= \ (QMiix) + Qau(-x)) , (17) 
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where Q Mi ^ } for i = 1,2,3, is the Gaussian state generated by the correlation matrix 
Mj(A), with A G [—1, l] 4 , and Mj(A) assumes one of the three forms: 

4 

M x (A) = 



M 2 (A) = 



M 3 (A) = 
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By construction, these states are convex-Gaussian, have null correlation matrix, and any 
convex combination of them has an expansion in terms of Majorana operators with non- 
vanishing coefficients only on the same correlators as Q a s{p)- The aim is to find for which 
p's it is possible to decompose g a ?,{p) as a convex sum of these type of convex- Gaussian 
states. 

In order to describe the convex hull of these states, first note that M2 and M3 are 
related to Mi by permutations (different choices of pairings): 



M 2 (A) 
M 3 (A) 



(P12 

(Pl3 



P l2 )M 1 {X) (P 12 ©P 12 ) 
P 13 )M X (A) (P 13 ©P 



where 



12 



/ 1 






o\ 




1 J 



p 



13 



/ 1 






13 J 




1 





\ 

1 







(18) 



Given that the transformation connecting the correlation matrices is formed by the 
direct sum of identical permutations, this transformation is in SO (8) - one always gets 
the signature of the permutation squared, and therefore the corresponding determinant 
is always 1. Rotations in SO of the correlation matrix induce unitary transformations 
on the level of the states [14] : Pi 2 t— > £/ 2 , P13 1— > U3. Bearing in mind this connection 
among the three classes of states above, we set up to determine the extreme points of 
the type 1 states. 

Lemma 3 (Extreme points for type 1 states). Let Si = conv{^i(A)|A G [—1, l] 4 }. Then 
the extreme points of Si are given by 



</>i(X) = o (\ x )( x \ + h x )h x \) > 



with x G {0, l} 4 and ->x is the bitwise negation of x. 
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Proof. It follows immediately from the expansion of the state Qi(X): 



£mi(a) 



2 t E II(w-ir fc A,o + n(i-(-ir*(-A,)) 

xe{0,l} 4 U=l k=l 

h E il(l-(-l) Xk h)(\x)(x\ + hx)hx\) 
xe{o,i} 4 k=i 

E ft (1- (-!)****) 



(19) 



x-e{o,i} 4 k=1 



□ 



Defining in a similar fashion £2, £3, 4>2(x) and ^3(0;), we can then define the convex 
hull of any convex combination of the three types of states as S = conv{iS>i, S 2 , S 3 }. 
More explicitly: 

37 37 



S = < ^^^(^(x) 



li(x) > 0,^^ 7i (a;) = lL (20) 

i=l x=0 1=1 x=0 ) 

where the summation in x is to be carried out over its binary expansion in four bits. 

The problem now reduces to determine for which values of p G [0, 1] the linear 
system Q a &{p) = Y^h=i J2l=o 1i( x )4 ) i(. x ), with {7j} a probability distribution, has a 
solution. Simple inspection shows that such Gaussian decomposition is always possible 
whenever p > 8/9, which is then an upper-bound on p*. 



5. Discussion: Classical simulation of fermionic quantum computation 

Contrary to bosonic linear optics, quantum computations based on fermionic linear 
optics (FLO) can be efficiently simulated by a classical computer even when augmented 
by measurements of fermionic number operators [T3| [T5l [T4"] . The simulation relies on 
the following facts: 

(i) efficient encoding of Gaussian states - a Gaussian state of 2m Majorana fermions 
is fully described by its correlation matrix M, with 0(m 2 ) elements. 

(ii) FLO transformations map Gaussian states onto Gaussian states, with efficient 
update rule - every FLO transformation V G £2™ induces a rotation R G 5*0 (2m) 
on the Majorana operators space, VciV^ = Yl^=i^-ij c j- This, in turn, induces the 
map M 1 — y RMR T on the 2m x 2m correlation matrix, and this update can be 
evaluated in 0(m 3 ) steps. 

(iii) efficient read out of measurement probability distributions - via Wick's theorem ([6]) 
the probability of measuring the population in fermionic modes can be done 
efficiently, as the Pfaffian of a 2m x 2m matrix can be evaluated in 0(m 3 ). 
Furthermore, number operator measurements project Gaussian states onto 
Gaussian states. 
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With these three ingredients, the evolution of any initial Gaussian state can 
be efficiently followed by a classical computer at any time in the computation, and 
the output distribution can be exactly evaluated. In references [HJ [16] the allowed 
transformations in (ii) were extended to noisy Gaussian maps. These are completely- 
positive channels, not necessarily FLO unitaries, that map Gaussian states onto 
Gaussian states, and as such still admit efficient classical simulation. 

Our results imply that if p > 8/9, Qasip) is convex- Gaussian and this allows 
for the classical simulation of the computation as follows. Given a noise strength 
p > 8/9, one finds the decomposition of Qasip) in S, by solving the linear system 
Qas{p) = Y^,=i YH x =o 7i( x )0j( a; ) ~~ this has to be done only once for a fixed p. As each 
4>i(x) is in itself a convex sum of two pure Gaussian states, this leads to a decomposition 
of g a s(p) in terms of at most 48 pure Gaussian states. Let pi be the probability 
associated with the z-th pure Gaussian state in such decomposition. Then, whenever 
a state Q a %{p) is requested in the computation, one samples from {p^} and chooses the 
corresponding pure Gaussian state. From this point onwards, the simulation follows the 
scheme summarized in (i), (ii) and (iii) above. 

The results in Fig. [1] leave some gaps in our understanding, in particular what is 
the precise value of p*. In order to improve the lower bounds on one could consider 
distillation processes beyond the restricted set of braiding operations, but that still map 
Gaussian states onto Gaussian states, leaving thus the set of convex-Gaussian states 
unchanged. The distillation protocol of depolarized GHZ^ states, which are isomorphic 
to Qasip), proposed by Diir and Cirac at first sight suggests that our classical simulation 
threshold is tight, as it can distill a perfect GHZ^ provided that p < 8/9. Nevertheless, 
their protocol uses non-FLO operations, as it assumes that one can locally create the 
GHZ4 and uses distilled 2-qubits maximally entangled pairs for distribution. The 
protocol devised by Murao et al. in [30], directly distills GHZ^ from depolarized copies 
of it for p < 0.7328. Despite the fact that this direct approach seems more related to our 
question, it employs gates that are not immediately translated into braiding operations 
or even to FLO transformations. Although the choice to distill via braiding operations 
is physically motivated - due to their topological protection-, a different distillation 
threshold from FLO operations would suggest that either Bravyi's distillation protocol 
can be improved to higher depolarizing noise rates, or that there exists a 'transition 
zone' in which the computation by braiding of Majorana fermions is neither quantumly- 
universal nor can be efficiently simulated classically. 
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Appendix A. Symmetric versus Gaussian-symmetric Subspace 

Inspired by [31], we define a 'FLO twirl' as the map S(g) = Iflo ®UgW ®W 

where j FLO dU is defined by taking the Haar measure over the real orthogonal matrices 
R induced by U. The map S(g) is normalized such that it is trace-preserving. 

Lemma 4. The projector onto the null-space of A = £\ q <g) q zs IIa = o = 
( 2 ™)iS>(|0, 0)(0, 0|) where |0) zs a (Gaussian) vacuum state with respect to some set of 
annihilation operators a^. Thus the states where ip is a pure fermionic Gaussian 

state span the null-space of A which implies that the null-space of A is a subspace of the 
symmetric subspace. 

Proof. In order to prove that Ha=o = ( 2 ™ ) «S(| 0, 0) (0, 0|), we note that both the l.h.s. 
and the r.h.s. are U ® [/-invariant where U is any FLO transformation. Thus instead 
of considering whether TrXILA=o = Tr(X( 2 ™)iS>(|0, 0)(0, 0|)) for any X, we can just 
consider the trace with respect to invariant objects S(X). It can be shown that A 1 for 
i = 0, 1, 2, ... , 2m (and linear combinations thereof) are the only invariants under the 
group U <8> U where U is FLO transformation [32]. Clearly TrA*( 2 ™)<S(|0, 0)(0, 0|) = 
= TrATlA = o for all i ^ while TrILA=o = ( 2 ™) fixes the overall prefactor. Having 
established the form of the projector, it follows directly that the states ip) for any 
Gaussian if) span the null-space (Assume this is false and hence a state in the null- 
space \x) = \Xin) + \Xout) where \x in ) is in the span of \ip,ip) while \\out) is w.l.o.g. 
orthogonal to any |V,V>- We have n A=0 |x) = |x> while (^)5(|0, 0)(0, 0|)| X ) = \Xin) 
arriving at a contradiction.) As ip) for Gaussian pure states if) span the null-space, 
and P\ip,tp) = \ip,ip) with P the SWAP operator, the null-space is a subspace of the 
symmetric subspace. □ 



Appendix B. SDP to detect non-convex- Gaussian states 

Recall the series of tests constructed to detect that if a state g e C 2m is not convex- 
Gaussian: 

Input: g G C<im and an integer n > 2 
Body: Is there g cxt E Cg£ 

S.t. P-KgextP-K = £?ext, VV G S n 
Tr2 ) ...,n^ext = Q 

Trfet = 1 
A k > l g ext = 0,Vk^l 

g cx t > o 

Output: yes (provide g ext ) or no 

Note that we have explicitly kept the permutation-symmetry constraint in this program. 

In this Appendix we convert this problem into a feasibility semi-definite 
programming problem. To do that we must show that the above problem can be cast 
as the standard form of SDP's, namely: 
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minimize c x 

subject to F + X %F% > 
Ax= b 

where x G R d is the unknown vector to be optimized over, c is a given vector in R d , 
and {Fi} i= Q^ ^cL are given symmetric matrices. For the equality constraint, A G W xd 
is a given matrix, with rank(v4) = p, and b G R p . The restriction on the rank of A 
demands that the linear system has at least one solution, and all the rows are linearly 
independent. 

The first step is to associate an Hermitian operator Mj with each term 
i k c ai c a2 . . . c a2k in the expression ([I]). The set {Mj} <j<22m-i-i spans all the even 
hermitian operators of Ci m - Then any state q G £2™ can be written as: 

2 2m-l_ 1 

g= a * M ^ (B- 1 ) 

k=0 

and we choose Mo = /. 

For definiteness, in what follows we restrict to the case n — 2. Extensions to larger 
n's follow the same structure and can be immediately constructed. That said, a general 
extension on C% m ® C% m can be written as: 

2 2m-l_ 1 

£? cx t = Pij M i® M 3i ( B -2) 

i,j=0 

with the coefficients G R, to be fixed by the constraints of our problem. The symmetry 
constraint can be taken directly into the parametrization of g ext by imposing: 

2 2m-l_ 1 

ft*t = Yl Pu M i® M i+ J2 P ij (M i ®M j + M j ®M i ). (B.3) 

i=0 0<i<j<2 2m-1 -l 

At this point, the initially 2 4m_2 coefficients of g ex t are reduced to 2 2m ~ 2 (l + 2 2m ~ 1 ). 
The imposition that Trx(f? e xt) = Q further demands: 

22m — 1 — 1 1 

j=0 j=0 

Therefore, (3 j = aj/2 m . In this way, further 2 2m_1 coefficients are determined, and thus 
remain 2 2m_3 (2 2m — 2) free parameters. Note that the normalization of g ext is already 
guaranteed as long as Tr(f?) = a = 1. 

The question whether there exists an assignment of the remaining free coefficients 
respecting the last two constraints can be immediately posed as a feasibility SDP: 

minimize 

subject to g ext > 

A&xt = 0. 
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(B.6) 



Forgetting for the moment that the correlators are in general complex, a direct 
correspondence with the SDP standard form can be done as follows: 

22m — 1 

F -»■ (3 0fi M ® M + ^ /3 j(M ® M,- + M,- <g) M ) (B.5) 

Mj®Mj, 1 < j < 2 2m - 1 - 1 
Mj®M k + M k ®Mj, 1 < j < k < 2 2m ~ 1 — 1 ' 

with x = ({/3jj}i<j<2 2m - 1 -i, {/3i,fc}i<j<fc<2 2m - 1 -i) T - 
The equality constraint then reads: 

2 2m_3 (2 2m -2) 

^ x t AF t = -AF . (B.7) 

i=i 

This can be cast in the form Ax = b, by writing A = ([Ai 71 j] 1<i <22m-3( 2 2m_2)), and 
b = — [AFo]. The notation [G] means a column matrix constructed by stacking each 
column of G below each other. 

To finish the construction we must take into account that we are possibly dealing 
with complex values. To do that we employ the following well known result. 



Lemma 5. Let Z = & e € dxd . Then Z>0 iffT(Z) 



Re(Z) -Im(Z) 
Im(Z) Re(Z) 



Proof. (=>) Let z = x + iy be an eigenvector of Z with eigenvalue A. Then it is easy to 
see that (x, y) T is an eigenvector of T(Z) with eigenvalue A. 

(<=) By the same token, let (x, y) be an eigenvector of T(Z) with eigenvalue A. 
Then z = x + iy is an eigenvector of Z with eigenvalue A. □ 

The final dimension of each F t is then 2 2m+1 x 2 2m+1 , and there are 2 2m " 3 (2 2m - 2) 
of such matrices. For the equality constraint we write 



Re(A) \ I Re{b) 
Im(A) j X ~ { Im(b) 



(B.8) 



The final matrix (Re(A), Im(A)) T has dimension 2 4m+1 x 2 2m ~ 3 (2 2m - 2). Most likely 
the rank of this matrix is not equal to the number of rows. Therefore, before plugging 
it into the SDP one must evaluate its reduced-row form (echelon form), and remove the 
all-zero lines. 

If one wishes to include the PPT test, as in Program[2l this can be done by rewriting 
the two positive constraints, g ext > and g^ t > 0, as a single one, g ext © > 0. One 
thus need to map each complex Fi in the positivity constraint to F[ = F{ © F^ 1 . The 
final dimension of each real F[ is then 2 2m+2 x 2 2m+2 . The equality constraint remains 
the same. 
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